1 DATA

1.1 Pop sive over time dataset

data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
  
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5
## 1          no
## 2          no
## 3          no
## 4          no
## 5          no
## 6          no
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571               Med                  5               6     8/25/21
## 572               Med                  5               6     8/25/21
## 573               Med                  5               6     8/25/21
## 574               Med                  5               6     8/25/21
## 575               Med                  5               6     8/25/21
## 576               Med                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5
## 571          no                     no         yes
## 572     extinct                    yes         yes
## 573          no                     no          no
## 574     extinct                    yes         yes
## 575          no                     no          no
## 576          no                     no          no
dim(data)
## [1] 576  23
summary(data)
##     Block            Old_Label            Label            Population       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  Genetic_Diversity  Generation_Parents Generation_Eggs Date_Census       
##  Length:576         Min.   :0.0        Min.   :1.0     Length:576        
##  Class :character   1st Qu.:1.0        1st Qu.:2.0     Class :character  
##  Mode  :character   Median :2.5        Median :3.5     Mode  :character  
##                     Mean   :2.5        Mean   :3.5                       
##                     3rd Qu.:4.0        3rd Qu.:5.0                       
##                     Max.   :5.0        Max.   :6.0                       
##                                                                          
##  Date_Start_Eggs    Date_End_Eggs       Image_Box1         Image_Box2       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     MC_Box1          MC_Box2       MC_Total_Adults     HC_Box1      
##  Min.   :  0.00   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
##  1st Qu.: 26.52   1st Qu.: 25.99   1st Qu.: 51.05   1st Qu.: 22.00  
##  Median : 58.87   Median : 54.29   Median :101.50   Median : 56.50  
##  Mean   : 75.88   Mean   : 70.57   Mean   :128.37   Mean   : 72.06  
##  3rd Qu.:106.93   3rd Qu.: 99.40   3rd Qu.:170.25   3rd Qu.:101.00  
##  Max.   :327.40   Max.   :299.00   Max.   :514.40   Max.   :288.00  
##  NA's   :142      NA's   :141      NA's   :5        NA's   :144     
##     HC_Box2       HC_Total_Adults   Nb_parents      Growth_rate    
##  Min.   :  0.00   Min.   :  0.0   Min.   :  0.00   Min.   :0.0000  
##  1st Qu.: 25.00   1st Qu.: 62.0   1st Qu.: 67.25   1st Qu.:0.4183  
##  Median : 52.00   Median :100.0   Median :100.00   Median :0.9386  
##  Mean   : 68.08   Mean   :132.1   Mean   :138.30   Mean   :1.4482  
##  3rd Qu.: 96.75   3rd Qu.:169.0   3rd Qu.:188.00   3rd Qu.:2.3587  
##  Max.   :290.00   Max.   :499.0   Max.   :499.00   Max.   :6.8571  
##  NA's   :142      NA's   :44      NA's   :122      NA's   :142     
##  First_Throw        Extinction_finalstatus Less_than_5       
##  Length:576         Length:576             Length:576        
##  Class :character   Class :character       Class :character  
##  Mode  :character   Mode  :character       Mode  :character  
##                                                              
##                                                              
##                                                              
## 
#Remove populations killed due to the pathogen
data_status  <- data.frame(data$Population,data$Extinction_finalstatus)
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)
## [1] 72 23
data<-data[data$Extinction_finalstatus!="kill",]

#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove 

#Check variables
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : chr  "Block3" "Block3" "Block3" "Block3" ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : chr  "High1" "High2" "High3" "High4" ...
##  $ Genetic_Diversity     : chr  "High" "High" "High" "High" ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: chr  "no" "no" "no" "no" ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)


#Order levels 
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)
## [1] "High"   "Medium" "Low"
data<-droplevels(data) #remove previous levels 
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
##  $ Genetic_Diversity     : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5 Nb_adults Lambda obs
## 1          no       100     NA   1
## 2          no       100     NA   2
## 3          no       100     NA   3
## 4          no       100     NA   4
## 5          no       100     NA   5
## 6          no       100     NA   6
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571            Medium                  5               6     8/25/21
## 572            Medium                  5               6     8/25/21
## 573            Medium                  5               6     8/25/21
## 574            Medium                  5               6     8/25/21
## 575            Medium                  5               6     8/25/21
## 576            Medium                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults   Lambda obs
## 571          no                     no         yes         4 0.500000 499
## 572     extinct                    yes         yes        NA       NA 500
## 573          no                     no          no        46 1.314286 501
## 574     extinct                    yes         yes        NA       NA 502
## 575          no                     no          no        56 1.400000 503
## 576          no                     no          no       133 3.243902 504
dim(data)
## [1] 504  26

1.2 Heterozygosity remain over time

#Upload He dataset
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]

#Merge dataset: add heterozygosity data to oridinal data 
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)

### Creation of 4 He: 
  # 1- He_start: He was remaining for each population when we started the experiment 
  # 2- He_lost_gen_t: He was lost only during this specfic generation during the experiment  (considering He=1 before this generation)
  # 3- He_remain: He was remaining after each generation (He_remain[Gen=6]=He_end)
  # 4- He_exp: He was lost during the entire experiment (considering He=1 at the beginning of the exp)
  # 5- He_end: He was remaining for each population when we finished the experiment (considering the real He at the beginning of the exp, He_start, not 1)


##### 
##### 1- He_start
#####
colnames(data)[27] <- "He_start"



##### 
##### 2- He_lost_gen_t
##### 
#Add He lost each generation
data$He_lost_gen_t <- 1 - (1/(2*data$Nb_adults))
#To consider extinct populations
is.na(data$He_lost_gen_t) <- sapply (data$He_lost_gen_t, is.infinite)


##### 
##### 3- He_remain
##### 
#Create new variable He remain per generation
data$He_remain <- NA
#Initiate with gen=1: He_start * He_lost_gen1
data$He_remain[data$Generation_Eggs=="1"] <- data$He_start[data$Generation_Eggs=="1"]*
  data$He_lost_gen_t[data$Generation_Eggs=="1"]

#Compute for all the other genrations, except the first one 
for(i in levels(data$Population)){
     ifelse(sum(data$Generation_Eggs[data$Population==i])!=21,print("ERROR"),print(i))
 for(j in 2:max(data$Generation_Eggs)){
   data$He_remain[data$Population==i&data$Generation_Eggs==j] <- 
     data$He_remain[data$Population==i&data$Generation_Eggs==j-1]*
     data$He_lost_gen_t[data$Population==i&data$Generation_Eggs==j]
   }
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
##### 
##### 4- He_exp
##### 
# He at the end of the experiment 
data$He_exp <- NA 

#Compute the HE at the end of the experiment 
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Eggs)!=21,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_adults))
#To consider extinct add this row: 
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
#Compute for the whole experiment 
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data$He_exp[data$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
##### 
##### 5- He_end
##### 
# Remaining He at the end of the experiment
data$He_end <- data$He_start * data$He_exp

#Save these values for Phenotyping 
data_he <- merge(x = data_he, 
                 y = data[data$Generation_Eggs==1,
                          c("Population","He_start",
                             "He_exp","He_end")], 
                 by = "Population", all.x=TRUE)

1.3 Phenotyping dataset

#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)


#Factors variables 
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <-   plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity, 
                                             levels = c("High", "Medium", "Low"))



#Add sum total 
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults  / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx

data_phenotyping_all <- data_phenotyping
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]




#Merge with He dataset
data_phenotyping <- merge(x = data_phenotyping, 
              y = data_he, by = "Population", all.x=TRUE)

#Clean if less than 30 parents 
pop_possible<-unique(data$Population[data$Generation_Eggs==5&data$Nb_adults>=30])
data_phenotyping <- data_phenotyping[data_phenotyping$Population %in% pop_possible,]

1.4 Old He remain

#Upload growth rate end of experiment
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]


### Additional: He remaining 
# New vector for the lost during exp
data_he$He_remain_exp <- NA

# He lost during exp
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))

#To consider extinct add this row: 
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)

temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
# Remaining He at the end of the experiment
data_he$He_end <- data_he$He_remain * data_he$He_remain_exp

# Remove kill population 
data_he <- data_he[!is.na(data_he$He_remain_exp),]
data_he <- droplevels(data_he)

## Merge dataset He
data_old_merged <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)

2 PLOT

2.1 Population size

PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Population,
                                colour = Genetic_Diversity)) + 
  geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
  geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

SUM_Popsize <- Rmisc::summarySE(data, 
                        measurevar = c("Nb_adults"),
                       groupvars = c("Generation_Eggs",
                                     "Genetic_Diversity"), 
                       na.rm=TRUE)
SUM_Popsize
##    Generation_Eggs Genetic_Diversity  N Nb_adults        sd        se       ci
## 1                1              High 27 100.00000   0.00000  0.000000  0.00000
## 2                1            Medium 23 100.00000   0.00000  0.000000  0.00000
## 3                1               Low 34 100.00000   0.00000  0.000000  0.00000
## 4                2              High 27 344.29630  73.77294 14.197610 29.18360
## 5                2            Medium 23 282.04348 100.75735 21.009360 43.57075
## 6                2               Low 34 282.61765  78.53006 13.467794 27.40043
## 7                3              High 27 151.85185  80.96899 15.582489 32.03027
## 8                3            Medium 23 118.73913  88.54491 18.462891 38.28969
## 9                3               Low 34 104.61765  85.13341 14.600260 29.70445
## 10               4              High 26 114.00000  80.20224 15.728954 32.39439
## 11               4            Medium 23  62.65217  64.50483 13.450188 27.89398
## 12               4               Low 34  46.38235  39.63241  6.796903 13.82840
## 13               5              High 27 103.62963  51.56486  9.923661 20.39838
## 14               5            Medium 20  61.50000  50.59904 11.314290 23.68108
## 15               5               Low 31  51.51613  46.13377  8.285870 16.92200
## 16               6              High 27 147.66667  42.60282  8.198916 16.85311
## 17               6            Medium 19  69.73684  46.02396 10.558620 22.18284
## 18               6               Low 28  77.03571  46.14428  8.720450 17.89289
# SUM_Popsize <- Rmisc::summarySE(data[data$Extinction_finalstatus=="no",],
#                         measurevar = c("Nb_adults"),
#                        groupvars = c("Generation_Eggs",
#                                      "Genetic_Diversity"),
#                        na.rm=TRUE)
# SUM_Popsize



PLOT_Pop_size_average <- PLOT_Pop_size + 
  geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 5, position = position_dodge(0.4)) +
  geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                ymin = Nb_adults-ci, ymax = Nb_adults+ci, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 1, width=.2,  position = position_dodge(0.4)) +
  geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
           linetype = "dashed", size = 1.5, position = position_dodge(0.4)) 
PLOT_Pop_size_average

# cowplot::save_plot(file = here::here("figures","PopulationSize_average.pdf"),   PLOT_Pop_size_average, base_height = 10/cm(1),           base_width = 15/cm(1), dpi = 610)






#Prediction with models 
# Test effect

data$gen_square <- data$Generation_Eggs^2
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], 
            family = "poisson")

#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(2,6, length = 100),each=3),
                       Block = "Block4",
                       gen_square = (rep(seq(2,6, length = 100),each=3))^2,
                       Estimates = NA)

filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)


#Predcit estimate minimun values 
tapply(filldata$Estimates, filldata$Genetic_Diversity, min)
##     High      Low   Medium 
## 85.85198 44.16308 51.99642
filldata[filldata$Estimates>=44.16&filldata$Estimates<=44.17,]
##     Genetic_Diversity Generation_Eggs  Block gen_square Estimates
## 198               Low        4.626263 Block4   21.40231  44.16308
filldata[filldata$Estimates>=51.99&filldata$Estimates<=52.00,]
##     Genetic_Diversity Generation_Eggs  Block gen_square Estimates
## 230            Medium        5.070707 Block4   25.71207  51.99642
filldata[filldata$Estimates>=85.85&filldata$Estimates<=85.86,]
##     Genetic_Diversity Generation_Eggs  Block gen_square Estimates
## 181              High        4.424242 Block4   19.57392  85.85198
## Change name vector 
vector_names <- c(`Low` = "Extreme-bottleneck",
                    `Medium` = "Severe-bottleneck",
                    `High` = "Diverse-source")


PLOT_Pop_size_predict<-ggplot2::ggplot(data[data$Extinction_finalstatus=="no",], 
                               aes(x = factor(Generation_Eggs), 
                                y = Nb_adults)) + 
  geom_point(aes(group = Population,
                 colour = Genetic_Diversity), 
             position = position_dodge(0.4), size = 0.7, alpha = 0.5) +
  geom_line(aes(group = Population,
                colour = Genetic_Diversity), 
            position = position_dodge(0.4), size = 0.05, alpha = 0.5) +
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.4) + 
  scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse-source","Severe-bottleneck","Extreme-bottleneck"), 
                    values = c("#00AFBB", "#E7B800", "#FC4E07")) +
  ylab("Population size") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size_predict

# 
# cowplot::save_plot(file = here::here("figures","PopulationSize_predict.pdf"),   PLOT_Pop_size_predict, base_height = 10/cm(1),           base_width = 15/cm(1), dpi = 610)

2.2 Extinction

#Percent of extinction
length(unique(data$Population[data$Genetic_Diversity=="Low"&
                                data$Extinction_finalstatus=="yes"]))/
  length(unique(data$Population[data$Genetic_Diversity=="Low"]))
## [1] 0.2352941
length(unique(data$Population[data$Genetic_Diversity=="Medium"&
                                data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="Medium"]))
## [1] 0.2173913
length(unique(data$Population[data$Genetic_Diversity=="High"&
                                data$Extinction_finalstatus=="yes"]))/
  length(unique(data$Population[data$Genetic_Diversity=="High"]))
## [1] 0
#Create data with percent
vector_names <- c(`Low` = "Extreme-bottleneck",
                    `Medium` = "Severe-bottleneck",
                    `High` = "Diverse-source")
f_labels = data.frame(Genetic_Diversity = c("High","Medium","Low"), 
                       label = c("Extinction = 0 %", "Extinction = 21.7 %", "Extinction = 23.5 %"))
f_labels$Genetic_Diversity <- factor(f_labels$Genetic_Diversity, 
                                     levels = c("High", "Medium", "Low"))






## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data,  aes(x = factor(Generation_Eggs), 
                                y = Nb_adults)) + 
  facet_wrap(~ Genetic_Diversity, ncol=3, labeller = ggplot2::as_labeller(vector_names)) + 
  #geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
  geom_line(aes(group = Population,
                colour = Extinction_finalstatus), 
            position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
  geom_text(x = 5, y = 420, 
            aes(label = label), 
            data = f_labels, 
            col="red", 
            size = 3, 
            vjust = 0) +
  scale_color_manual(values = c("black", "red")) +
  ylab("Population size") +
  xlab("Generation") +
  theme(legend.position = "none") +
  theme_LO
PLOT_Extinction

cowplot::save_plot(here::here("figures","Extinction.pdf"),   PLOT_Extinction, base_height = 8/cm(1),
          base_width = 20/cm(1), dpi = 610)
# 




PLOT_Extinction 

2.3 Growth rate G2

tapply(data[data$Generation_Eggs==2,]$Lambda,data[data$Generation_Eggs==2,]$Genetic_Diversity,mean)
##     High   Medium      Low 
## 3.442963 2.820435 2.826176
PLOT_Growth<-ggplot2::ggplot(data[data$Generation_Eggs==2,], aes(x = Genetic_Diversity, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
  geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Genetic diversity") +
  theme_LO
PLOT_Growth

# 
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"),   PLOT_Growth, base_height = 10/cm(1),
#           base_width = 8/cm(1), dpi = 610)

2.4 Life stage proportion

SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_adults"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_adults
##   Genetic_Diversity   Week   N  p_adults         sd          se         ci
## 1              High 4-week  50 0.6508098 0.17776618 0.025139934 0.05052059
## 2              High 5-week 105 0.9632040 0.05259208 0.005132461 0.01017786
## 3            Medium 4-week  24 0.7121270 0.15722597 0.032093616 0.06639070
## 4            Medium 5-week  49 0.9563136 0.06160975 0.008801393 0.01769639
## 5               Low 4-week  30 0.6568606 0.19029917 0.034743716 0.07105888
## 6               Low 5-week  59 0.9582021 0.07845154 0.010213520 0.02044458
SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_pupae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_pupae
##   Genetic_Diversity   Week   N    p_pupae         sd          se          ci
## 1              High 4-week  50 0.28825555 0.15629878 0.022103985 0.044419621
## 2              High 5-week 105 0.01458769 0.02407445 0.002349426 0.004658999
## 3            Medium 4-week  24 0.23019148 0.14506173 0.029610601 0.061254195
## 4            Medium 5-week  49 0.02283167 0.03485130 0.004978757 0.010010462
## 5               Low 4-week  30 0.29380227 0.17562943 0.032065401 0.065581108
## 6               Low 5-week  59 0.01675308 0.02459640 0.003202178 0.006409856
SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_larvae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_larvae
##   Genetic_Diversity   Week   N   p_larvae         sd          se          ci
## 1              High 4-week  50 0.06093466 0.04274850 0.006045551 0.012148989
## 2              High 5-week 105 0.02220828 0.04324043 0.004219833 0.008368088
## 3            Medium 4-week  24 0.05768150 0.05999020 0.012245449 0.025331640
## 4            Medium 5-week  49 0.02085474 0.04462923 0.006375605 0.012819012
## 5               Low 4-week  30 0.04933713 0.08198954 0.014969174 0.030615398
## 6               Low 5-week  59 0.02504479 0.06691319 0.008711355 0.017437671
SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
                       p_pupae=SUM_Prop_pupae$p_pupae,
                       p_larvae=SUM_Prop_larvae$p_larvae)


SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop
##    Genetic_Diversity   Week   N    Stage Proportion
## 1               High 4-week  50 p_adults 0.65080979
## 2               High 5-week 105 p_adults 0.96320403
## 3             Medium 4-week  24 p_adults 0.71212702
## 4             Medium 5-week  49 p_adults 0.95631359
## 5                Low 4-week  30 p_adults 0.65686059
## 6                Low 5-week  59 p_adults 0.95820213
## 7               High 4-week  50  p_pupae 0.28825555
## 8               High 5-week 105  p_pupae 0.01458769
## 9             Medium 4-week  24  p_pupae 0.23019148
## 10            Medium 5-week  49  p_pupae 0.02283167
## 11               Low 4-week  30  p_pupae 0.29380227
## 12               Low 5-week  59  p_pupae 0.01675308
## 13              High 4-week  50 p_larvae 0.06093466
## 14              High 5-week 105 p_larvae 0.02220828
## 15            Medium 4-week  24 p_larvae 0.05768150
## 16            Medium 5-week  49 p_larvae 0.02085474
## 17               Low 4-week  30 p_larvae 0.04933713
## 18               Low 5-week  59 p_larvae 0.02504479
SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))





PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity, 
                                                  y = Proportion, fill = Stage)) +
  facet_wrap(~ Week, ncol=2) +
  geom_bar(stat="identity") + 
  xlab("Genetic history") +
  labs(color="Stage of individuals") +
  ylab("Proportion of each stage") +
  scale_fill_manual(values= c("#619CFF","#F8766D","#00BA38"), 
                    breaks = c("p_adults", "p_pupae", "p_larvae"),
                    labels = c( "Adult", "Pupae","Larvae")) + 
  ggplot2::scale_x_discrete(breaks=c("High", "Medium", "Low"), 
                            labels=c("Diverse-\nsource", 
                                     "Severe-\nbottleneck", 
                                      "Extreme-\nbottleneck")) + 
  theme_LO
PLOT_prop

# cowplot::save_plot(here::here("figures","Life_stage_proportion.pdf"),   PLOT_prop,
#           base_height = 10/cm(1), base_width = 18/cm(1), dpi = 610)

2.5 Remaining heterozygosity

# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_phenotyping,
                            measurevar = c("Lambda"),
                            groupvars = c("Genetic_Diversity",
                                          "He_end",
                                          "Population"), 
                            na.rm=TRUE)


PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_end,
                                                      y = Lambda, 
                                                      colour = Genetic_Diversity, 
                                                      size = N)) + 
    geom_point(alpha = 0.8) +
    ylab("Growth rate") +
    xlab("Expected heterozygosity") +
    scale_color_manual(name="Genetic history",
                     breaks = c("High", "Medium", "Low"), 
                    labels = c("Diverse-source","Severe-bottleneck","Extreme-bottleneck"), 
                    values = c("#00AFBB", "#E7B800", "#FC4E07")) + 
    labs(size="Replicates") + 
    theme_LO 
PLOT_He_Final

# 
#    
# cowplot::save_plot(here::here("figures","Heterozygosity_Growthrate.pdf"),
#                    PLOT_He_Final,
#           base_height = 12/cm(1), base_width = 16/cm(1), dpi = 610)
# 




# # ALL WITHOUT PSEUDO
# SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
#                             measurevar = c("Lambda"),
#                        groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_end",
#                                      "Population"), 
#                        na.rm=TRUE)
# 
# SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_end)
# 
# SUM_POP_He_ALL <- SUM_POP_He_ALL[!is.na(SUM_POP_He_ALL$He_end),]
# SUM_POP_He_ALL <- SUM_POP_He_ALL[SUM_POP_He_ALL$HE!="-Inf",]
# 
# 
# facet_names <- c('Initial'="Before rescue experiment",'Final'="After rescue experiment")
# 
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
#                          colour = Genetic_Diversity, 
#                          shape = Phenotyping)) +
#   #facet_wrap(~ Phenotyping, ncol=1) +
#   geom_point(size = 3, alpha = 0.8) +
#   scale_shape_manual(values = c(1, 19)) + 
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO 
# PLOT_He_ALL
# 
# 
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
#                          colour = Genetic_Diversity, 
#                          shape = Phenotyping)) +
#   facet_wrap(~ Phenotyping, ncol=1) +
#   geom_point(size = 3, alpha = 0.8) +
#   scale_shape_manual(values = c(1, 19)) + 
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO
# PLOT_He_ALL
# 
# 
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
#                          colour = Genetic_Diversity)) +
#   facet_wrap(~ Phenotyping, ncol=1) +
#   geom_point(size = 3, alpha = 0.8) +
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO
# PLOT_He_ALL

2.6 He over time

#Compute for all the other genrations, except the first one 
data$ID_extinction <- "extant"

for(i in unique(data$Population[data$Extinction_finalstatus=="yes"])){
    gen_all <- data$Generation_Eggs[data$Population==i&!is.na(data$He_remain)]
    maxgen <- max(gen_all)
    data$ID_extinction[data$Population==i&data$Generation_Eggs==maxgen] <- "willextinct" 
}

data[data$ID_extinction=="willextinct",]
##     Population  Block         Old_Label                  Label
## 205      Low16 Block4             B4-O1 7/14 BE B4 | G5 Low 16
## 273      Low27 Block5             B5-B1 6/16 BE B5 | G4 Low 27
## 278      Low28 Block5             B5_E1 5/12 BE B5 | G3 Low 28
## 299      Low30 Block5             B5-D1 6/16 BE B5 | G4 Low 30
## 303      Low31 Block5         B(2)5-P1  6/16 BE B5 | G4 Low 31
## 310      Low32 Block5         B(2)5_Q1  5/12 BE B5 | G3 Low 32
## 317      Low33 Block5         B(2)5-T1  7/21 BE B5 | G5 Low 33
## 336      Low36 Block5         B(2)5-X1  6/16 BE B5 | G4 Low 36
## 402      Med14 Block4   B4-Backup Fam L  6/9 BE B4 | G4 Med 14
## 409      Med17 Block5   B5_Backup_Fam_F 5/12 BE B5 | G3 Med 17
## 437      Med20 Block5 B5 - Backup Fam I 6/16 BE B5 | G4 Med 20
## 449      Med22 Block5   B5_Backup_Fam_B 5/12 BE B5 | G3 Med 22
## 479       Med5 Block3 B3 - Backup Fam E   7/7 BE B3 | G5 Med 5
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 205               Low                  4               5     7/14/21
## 273               Low                  3               4     6/16/21
## 278               Low                  2               3     5/12/21
## 299               Low                  3               4     6/16/21
## 303               Low                  3               4     6/16/21
## 310               Low                  2               3     5/12/21
## 317               Low                  4               5     7/21/21
## 336               Low                  3               4     6/16/21
## 402            Medium                  3               4      6/9/21
## 409            Medium                  2               3     5/12/21
## 437            Medium                  3               4     6/16/21
## 449            Medium                  2               3     5/12/21
## 479            Medium                  4               5      7/7/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2   MC_Box1   MC_Box2
## 205         7/14/21       7/15/21       <NA>   DSC_0994  0.000000  2.000000
## 273         6/16/21       6/17/21   DSC_0526   DSC_0527  5.090878  0.000000
## 278         5/12/21       5/13/21   DSC_0918   DSC_0919  3.085768  3.638847
## 299         6/16/21       6/17/21   DSC_0559       <NA>  1.383920  0.000000
## 303         6/16/21       6/17/21       <NA>   DSC_0541  0.000000  1.000000
## 310         5/12/21       5/13/21   DSC_0953   DSC_0954 88.254170 70.952124
## 317         7/21/21       7/22/21   DSC_0117   DSC_0118  1.000000  1.000000
## 336         6/16/21       6/17/21   DSC_0540       <NA>  1.000000  0.000000
## 402          6/9/21       6/10/21   DSC_0376   DSC_0377 10.974990  1.000000
## 409         5/12/21       5/13/21   DSC_0916   DSC_0917  1.543122  1.084327
## 437         6/16/21       6/17/21   DSC_0542       <NA>  1.000000  0.000000
## 449         5/12/21       5/13/21   DSC_0922   DSC_0923  9.168447  6.251069
## 479          7/7/21        7/8/21       <NA>   DSC_0830        NA        NA
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 205             2.0       0       2               2          5 0.400000000
## 273             5.1       3       0               3         19 0.157894737
## 278             6.7       1       1               2        374 0.005347594
## 299             1.4       1       0               1        146 0.006849315
## 303             1.0       1       1               2        272 0.007352941
## 310           159.2      71      56             127        276 0.460144928
## 317             2.0       1       1               2          3 0.666666667
## 336             1.0       1       0               1         22 0.045454545
## 402            12.0       7       1               8        138 0.057971014
## 409             2.6       3       2               5        416 0.012019231
## 437             1.0       1       0               1         20 0.050000000
## 449            15.4      10       6              16        200 0.080000000
## 479             0.0      NA      NA               1         11 0.090909091
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults      Lambda obs
## 205          no                    yes         yes         2 0.400000000 378
## 273          no                    yes         yes         3 0.157894737 319
## 278          no                    yes         yes         2 0.005347594 236
## 299          no                    yes         yes         1 0.006849315 322
## 303          no                    yes         yes         2 0.007352941 323
## 310          no                    yes         yes       127 0.460144928 240
## 317          no                    yes         yes         2 0.666666667 409
## 336          no                    yes         yes         1 0.045454545 328
## 402          no                    yes         yes         8 0.057971014 308
## 409          no                    yes         yes         5 0.012019231 245
## 437          no                    yes         yes         1 0.050000000 332
## 449          no                    yes         yes        16 0.080000000 250
## 479          no                    yes         yes         1 0.090909091 359
##      He_start He_lost_gen_t He_remain    He_exp    He_end gen_square
## 205 0.5544299     0.7500000 0.3575672 0.6449276 0.3575672         25
## 273 0.5367998     0.8333333 0.4324691 0.8056432 0.4324691         16
## 278 0.5386585     0.7500000 0.4014365 0.7452523 0.4014365          9
## 299 0.5540444     0.5000000 0.2742819 0.4950540 0.2742819         16
## 303 0.5532775     0.7500000 0.4116634 0.7440450 0.4116634         16
## 310 0.5528880     0.9960630 0.5469651 0.9892872 0.5469651          9
## 317 0.5523333     0.7500000 0.3359893 0.6083089 0.3359893         25
## 336 0.5427371     0.5000000 0.2635269 0.4855518 0.2635269         16
## 402 0.6802331     0.9375000 0.6315974 0.9285014 0.6315974         16
## 409 0.6825509     0.9000000 0.6104897 0.8944237 0.6104897          9
## 437 0.6819650     0.5000000 0.3302071 0.4841994 0.3302071         16
## 449 0.6809168     0.9687500 0.6546991 0.9614965 0.6546991          9
## 479 0.6807065     0.5000000 0.3209456 0.4714889 0.3209456         25
##     ID_extinction
## 205   willextinct
## 273   willextinct
## 278   willextinct
## 299   willextinct
## 303   willextinct
## 310   willextinct
## 317   willextinct
## 336   willextinct
## 402   willextinct
## 409   willextinct
## 437   willextinct
## 449   willextinct
## 479   willextinct
dim(data[data$ID_extinction=="willextinct",])
## [1] 13 33
data[data$ID_extinction=="willextinct","Population"]
##  [1] Low16 Low27 Low28 Low30 Low31 Low32 Low33 Low36 Med14 Med17 Med20 Med22
## [13] Med5 
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9
data$Population[data$Genetic_Diversity=="Medium"&data$Extinction_finalstatus=="yes"]
##  [1] Med14 Med14 Med14 Med14 Med14 Med14 Med17 Med17 Med17 Med17 Med17 Med17
## [13] Med20 Med20 Med20 Med20 Med20 Med20 Med22 Med22 Med22 Med22 Med22 Med22
## [25] Med5  Med5  Med5  Med5  Med5  Med5 
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9
PLOT_He_extinction <-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = He_remain, 
                                group = Population,
                                colour = Extinction_finalstatus, 
                                shape = ID_extinction)) + 
  geom_point(position = position_dodge(0.5), size = 3, alpha = 0.6) +
  geom_line(position = position_dodge(0.5), size = 0.25, alpha = 0.4) +
  ylab("Expected heterozygozity") +
  scale_color_manual(values = c("black", "red")) +
  scale_shape_manual(values = c(16, 13), guide=FALSE) +
  labs(color="Extinct populations") +
  xlab("Generation") +
  theme_LO
PLOT_He_extinction

# 
# cowplot::save_plot(here::here("figures","Heterozygosity_over_time.pdf"),
#                    PLOT_He_extinction,
#           base_height = 12/cm(1), base_width = 14/cm(1), dpi = 610)

3 ANALYSIS

3.1 Probability of extinction

############ 
###### Clean dataset
############ 
#Prepare clean dataset
data_proba_extinction <- data[data$Generation_Eggs==6,c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)



############ 
###### Analysis
############ 
#Analysis
mod0 <- glm(y ~ Genetic_Diversity + Block,  
      data = data_proba_extinction, family = binomial)


summary(mod0)
## 
## Call:
## glm(formula = y ~ Genetic_Diversity + Block, family = binomial, 
##     data = data_proba_extinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.26011  -0.44258  -0.30957  -0.00005   2.47474  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -21.3144  1855.7435  -0.011  0.99084   
## Genetic_DiversityMedium   18.3002  1855.7433   0.010  0.99213   
## Genetic_DiversityLow      18.5062  1855.7432   0.010  0.99204   
## BlockBlock4                0.7402     1.2714   0.582  0.56046   
## BlockBlock5                3.0006     1.1265   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 72.388  on 83  degrees of freedom
## Residual deviance: 46.833  on 79  degrees of freedom
## AIC: 56.833
## 
## Number of Fisher Scoring iterations: 18
#Test genetic diversity effect
mod1 <- glm(y ~ Block ,  
      data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1        81     58.903                        
## 2        79     46.833  2   12.069 0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#We keep genetic diversity 

# Perform the lrtest 
(A <- logLik(mod1))
## 'log Lik.' -29.45146 (df=3)
(B <- logLik(mod0))
## 'log Lik.' -23.41669 (df=5)
#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))
## [1] 12.06954
#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))
## [1] 0.002394043
lmtest::lrtest(mod1, mod0)
## Likelihood ratio test
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   #Df  LogLik Df  Chisq Pr(>Chisq)   
## 1   3 -29.451                        
## 2   5 -23.417  2 12.069   0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,  
      data = data_proba_extinction, family = binomial)
summary(mod0)
## 
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial, 
##     data = data_proba_extinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.26011  -0.44258  -0.30957  -0.00005   2.47474  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## Genetic_DiversityHigh    -21.3144  1855.7435  -0.011  0.99084   
## Genetic_DiversityMedium   -3.0142     1.1293  -2.669  0.00760 **
## Genetic_DiversityLow      -2.8082     1.0659  -2.635  0.00843 **
## BlockBlock4                0.7402     1.2714   0.582  0.56046   
## BlockBlock5                3.0006     1.1265   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.449  on 84  degrees of freedom
## Residual deviance:  46.833  on 79  degrees of freedom
## AIC: 56.833
## 
## Number of Fisher Scoring iterations: 18
############ 
###### Predicted value
############ 
#Extract probability of extinction 
(p_high <- plogis(-21.3144))
## [1] 5.536989e-10
(p_med <- plogis(-3.0142))
## [1] 0.04678847
(p_low <- plogis(-2.8082))
## [1] 0.05688267
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
                                   Block = "Block4")

data_predict_extinct$predict <- predict(mod0, type = "response",
        re.form = NA,
        newdata = data_predict_extinct)
        
data_predict_extinct
##   Genetic_Diversity  Block      predict
## 1              High Block4 1.160723e-09
## 2            Medium Block4 9.329490e-02
## 3               Low Block4 1.122446e-01
p_high
## [1] 5.536989e-10
p_med
## [1] 0.04678847
p_low
## [1] 0.05688267
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE  df asymp.LCL asymp.UCL
##  High              -20.07 1855.743 Inf  -4451.10  4410.961
##  Medium             -1.77    0.650 Inf     -3.32    -0.216
##  Low                -1.56    0.532 Inf     -2.83    -0.290
## 
## Results are averaged over the levels of: Block 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate       SE  df z.ratio p.value
##  High - Medium  -18.300 1855.743 Inf -0.010  0.9999 
##  High - Low     -18.506 1855.743 Inf -0.010  0.9999 
##  Medium - Low    -0.206    0.751 Inf -0.274  0.9594 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.2 Population size G6

############ 
###### Analysis
############ 
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),  
#       data = data[data$Generation_Eggs==6,], family = "poisson")


mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])

mod1 <- lm(log(Nb_adults) ~ Block , 
             data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])

anova(mod0, mod1) # Effect of genetic diversity 
## Analysis of Variance Table
## 
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     66 29.706                                  
## 2     68 42.118 -2   -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])


summary(mod0)
## 
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data[data$Generation_Eggs == 
##     6 & data$Extinction_finalstatus == "no", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2417 -0.2947  0.1088  0.4381  1.2624 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              5.11532    0.17980  28.451  < 2e-16 ***
## Genetic_DiversityMedium -0.94813    0.20540  -4.616 1.86e-05 ***
## Genetic_DiversityLow    -0.80532    0.18719  -4.302 5.72e-05 ***
## BlockBlock4             -0.02077    0.18447  -0.113   0.9107    
## BlockBlock5             -0.53922    0.21414  -2.518   0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared:  0.3362, Adjusted R-squared:  0.2959 
## F-statistic: 8.356 on 4 and 66 DF,  p-value: 1.623e-05
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE df lower.CL upper.CL
##  High                4.93 0.130 66     4.61     5.25
##  Medium              3.98 0.159 66     3.59     4.37
##  Low                 4.12 0.136 66     3.79     4.46
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE df t.ratio p.value
##  High - Medium    0.948 0.205 66  4.616  0.0001 
##  High - Low       0.805 0.187 66  4.302  0.0002 
##  Medium - Low    -0.143 0.207 66 -0.690  0.7704 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.3 Population size over time

####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY WITHOUT GEN 1  
####################### 
###################### REMOVE GEN 1 ##########################33

#If we dont consider Gen=1 
PLOT_Pop_size_average

fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2,])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), data = data[data$Generation_Eggs>=2,])

fit <- glm(Nb_adults ~ Generation_Eggs, 
           data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit6 <- glm(Nb_adults~exp(Generation_Eggs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit7 <- glm(Nb_adults~log(Generation_Eggs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")





#generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(1,6, length=100)
# plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
#      log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
# lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
# lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
# lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
# lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
# lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
# lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
     data[data$Generation_Eggs>=2,]$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
rsq::rsq(fit,adj=TRUE)
## [1] 0.4442203
rsq::rsq(fit2,adj=TRUE)
## [1] 0.5938535
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5928435
rsq::rsq(fit4,adj=TRUE)
## [1] 0.592052
rsq::rsq(fit5,adj=TRUE)
## [1] 0.592052
rsq::rsq(fit6,adj=TRUE)
## [1] 0.121301
rsq::rsq(fit7,adj=TRUE)
## [1] 0.5205104
# Best model is fit 2


data$gen_square <- data$Generation_Eggs^2
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~Generation_Eggs + gen_square, 
            data = data[data$Generation_Eggs>=2,], family = "poisson")

# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
#Add the interaction with Genetic diversity
fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
             gen_square*Genetic_Diversity,
           data = data, family = "poisson")


summary(fit2)
## 
## Call:
## glm(formula = Nb_adults ~ Generation_Eggs * Genetic_Diversity + 
##     gen_square * Genetic_Diversity, family = "poisson", data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -15.961   -6.480   -2.662    5.254   21.115  
## 
## Coefficients:
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              5.070332   0.026463 191.604  < 2e-16
## Generation_Eggs                          0.156559   0.017811   8.790  < 2e-16
## Genetic_DiversityMedium                 -0.164697   0.042046  -3.917 8.96e-05
## Genetic_DiversityLow                     0.024718   0.037632   0.657 0.511290
## gen_square                              -0.036751   0.002566 -14.324  < 2e-16
## Generation_Eggs:Genetic_DiversityMedium  0.080423   0.029501   2.726 0.006409
## Generation_Eggs:Genetic_DiversityLow    -0.096063   0.026425  -3.635 0.000278
## Genetic_DiversityMedium:gen_square      -0.035166   0.004448  -7.905 2.67e-15
## Genetic_DiversityLow:gen_square         -0.009935   0.003956  -2.511 0.012038
##                                            
## (Intercept)                             ***
## Generation_Eggs                         ***
## Genetic_DiversityMedium                 ***
## Genetic_DiversityLow                       
## gen_square                              ***
## Generation_Eggs:Genetic_DiversityMedium ** 
## Generation_Eggs:Genetic_DiversityLow    ***
## Genetic_DiversityMedium:gen_square      ***
## Genetic_DiversityLow:gen_square         *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 39168  on 486  degrees of freedom
## Residual deviance: 30719  on 478  degrees of freedom
##   (17 observations deleted due to missingness)
## AIC: 33747
## 
## Number of Fisher Scoring iterations: 5
####### Add genetic diversity 

#Test oversdispersion
mod0 <- glm(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
              Block , 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
sqrt(deviance(mod0)/(nobs(mod0)-length(coef(mod0))))
## [1] 6.395475
sigma(mod0)
## [1] 6.395475
# ccl: overdispersion


mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
blmeco::dispersion_glmer(mod1)
## [1] 1.076696
sigma(mod1)
## [1] 1
#Convergence issue 
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.004551091
# Test effect
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
anova(mod1, mod2, test ="Chisq")  
## Data: data[data$Generation_Eggs >= 2, ]
## Models:
## mod2: Nb_adults ~ Generation_Eggs + gen_square + Genetic_Diversity + 
## mod2:     Block + (1 | obs)
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## mod1:     Genetic_Diversity + Block + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## mod2    8 4669.1 4701.1 -2326.6   4653.1                        
## mod1   12 4660.8 4708.8 -2318.4   4636.8 16.376  4   0.002554 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(2,6, length = 100),each=3),
                       Block = "Block4",
                       gen_square = (rep(seq(2,6, length = 100),each=3))^2,
                       Estimates = NA)

filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data[data$Generation_Eggs>=2,]) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

emmeans::emmeans(mod1, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean     SE  df asymp.LCL asymp.UCL
##  High                4.88 0.0790 Inf      4.69      5.07
##  Medium              4.19 0.0885 Inf      3.98      4.40
##  Low                 4.04 0.0733 Inf      3.87      4.22
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.688 0.119 Inf 5.786   <.0001 
##  High - Low       0.835 0.107 Inf 7.778   <.0001 
##  Medium - Low     0.147 0.115 Inf 1.287   0.4028 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
###################################################################
###################################################################
###################################################################
####################  INCLUDING GEN 1 ###########################
PLOT_Pop_size_average

fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")


#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
     data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
summary(fit)$adj.r.squared
## NULL
summary(fit2)$adj.r.squared
## NULL
summary(fit3)$adj.r.squared
## NULL
summary(fit4)$adj.r.squared # best r square
## NULL
summary(fit5)$adj.r.squared 
## NULL
summary(fit6)$adj.r.squared 
## NULL
summary(fit7)$adj.r.squared 
## NULL
rsq::rsq(fit,adj=TRUE)
## [1] 0.1067643
rsq::rsq(fit2,adj=TRUE)
## [1] 0.1486733
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5240033
rsq::rsq(fit4,adj=TRUE)
## [1] 0.5980368
rsq::rsq(fit5,adj=TRUE)
## [1] 0.5989426
rsq::rsq(fit6,adj=TRUE)
## [1] 0.07185974
rsq::rsq(fit7,adj=TRUE)
## [1] 0.04937507
#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")



####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY and GEN 1 
####################### 

#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)
## [1] 5.760264
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)
## [1] 1.083377
#Convergence issue 
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.06384154
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(1,6, length = 100),each=3),
                       Block = "Block4", 
                       Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity + 
#              gen_square*Genetic_Diversity, 
#            data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity + 
#              gen_square, 
#            data = data[data$Generation_Eggs>=2,])

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## mod1    9 5611.2 5648.9 -2796.6   5593.2                        
## mod    17 5594.8 5666.0 -2780.4   5560.8 32.47  8  7.672e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Signifcant interaction 
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5594.8   5666.0  -2780.4   5560.8      470 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.03685 -0.03969  0.00800  0.05695  0.26801 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.64646  0.8040  
##  Block  (Intercept) 0.07945  0.2819  
## Number of obs: 487, groups:  obs, 487; Block, 3
## 
## Fixed effects:
##                                                                       Estimate
## (Intercept)                                                          -1.993477
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.858800
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.147603
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.942314
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.058998
## Genetic_DiversityMedium                                              -0.352704
## Genetic_DiversityLow                                                 -0.787854
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  0.772686
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.503761
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.095237
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.006018
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     1.489498
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    -0.833931
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.138807
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    -0.007028
##                                                                      Std. Error
## (Intercept)                                                            1.238244
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           1.985855
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                           1.032426
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           0.213266
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                           0.015171
## Genetic_DiversityMedium                                                1.833248
## Genetic_DiversityLow                                                   1.630486
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   2.976126
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   1.553127
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.322073
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.022992
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      2.637956
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      1.373155
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.284239
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.020268
##                                                                      z value
## (Intercept)                                                           -1.610
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           5.468
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -4.986
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           4.418
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -3.889
## Genetic_DiversityMedium                                               -0.192
## Genetic_DiversityLow                                                  -0.483
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.260
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  -0.324
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.296
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  -0.262
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.565
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     -0.607
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.488
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     -0.347
##                                                                      Pr(>|z|)
## (Intercept)                                                          0.107415
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         4.55e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         6.17e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         9.94e-06
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         0.000101
## Genetic_DiversityMedium                                              0.847434
## Genetic_DiversityLow                                                 0.628952
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.795151
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.745671
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.767459
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.793538
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    0.572318
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    0.543645
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    0.625306
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    0.728786
##                                                                         
## (Intercept)                                                             
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.398806 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
summary(mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5611.2   5648.9  -2796.6   5593.2      478 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.05397 -0.05185  0.01426  0.07054  0.21832 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.69318  0.8326  
##  Block  (Intercept) 0.07665  0.2769  
## Number of obs: 487, groups:  obs, 487; Block, 3
## 
## Fixed effects:
##                                              Estimate Std. Error z value
## (Intercept)                                  -2.01996    0.81642  -2.474
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.74179    1.30992   8.964
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.66290    0.68789  -8.232
## stats::poly(Generation_Eggs, 4, raw = TRUE)3  1.03332    0.14306   7.223
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.06412    0.01022  -6.274
## Genetic_DiversityMedium                      -0.55984    0.10012  -5.592
## Genetic_DiversityLow                         -0.67366    0.09048  -7.445
##                                              Pr(>|z|)    
## (Intercept)                                    0.0134 *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.08e-13 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 3.52e-10 ***
## Genetic_DiversityMedium                      2.25e-08 ***
## Genetic_DiversityLow                         9.69e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                    (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.965                                      
## s::(G_E,4,r=TRUE)2  0.941 -0.993                               
## s::(G_E,4,r=TRUE)3 -0.915  0.976             -0.995            
## s::(G_E,4,r=TRUE)4  0.889 -0.956              0.984            
## Gntc_DvrstM        -0.062  0.007             -0.008            
## Gntc_DvrstL        -0.065  0.004             -0.005            
##                    s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1                                             
## s::(G_E,4,r=TRUE)2                                             
## s::(G_E,4,r=TRUE)3                                             
## s::(G_E,4,r=TRUE)4 -0.997                                      
## Gntc_DvrstM         0.008             -0.009                   
## Gntc_DvrstL         0.005             -0.006              0.494
## convergence code: 0
## Model failed to converge with max|grad| = 0.053402 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.55 0.206 Inf      4.06      5.04
##  Medium              3.93 0.216 Inf      3.42      4.45
##  Low                 3.69 0.201 Inf      3.21      4.17
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.619 0.187 Inf 3.310   0.0027 
##  High - Low       0.861 0.168 Inf 5.115   <.0001 
##  Medium - Low     0.242 0.184 Inf 1.319   0.3844 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
####################### 
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
####################### 
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
            data = data[data$Extinction_finalstatus=="no",], family = "poisson")

filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")

PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates_Withoutextpop, 
                               colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
  geom_point(data = data[data$Extinction_finalstatus=="no",], 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), 
             data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), 
              data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod1    9 4759.9 4796.3 -2370.9   4741.9                       
## mod    17 4756.3 4825.1 -2361.1   4722.3 19.599  8    0.01196 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data[data$Extinction_finalstatus == "no", ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   4756.3   4825.1  -2361.1   4722.3      407 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57063 -0.04879  0.00393  0.06967  0.29650 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.36624  0.6052  
##  Block  (Intercept) 0.02456  0.1567  
## Number of obs: 424, groups:  obs, 424; Block, 3
## 
## Fixed effects:
##                                                                        Estimate
## (Intercept)                                                          -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.0583188
## Genetic_DiversityMedium                                               1.1892488
## Genetic_DiversityLow                                                  0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0006044
##                                                                      Std. Error
## (Intercept)                                                           0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                          1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          0.0120893
## Genetic_DiversityMedium                                               1.5049889
## Genetic_DiversityLow                                                  1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0173582
##                                                                      z value
## (Intercept)                                                           -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -4.824
## Genetic_DiversityMedium                                                0.790
## Genetic_DiversityLow                                                   0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.035
##                                                                      Pr(>|z|)
## (Intercept)                                                            0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         1.41e-06
## Genetic_DiversityMedium                                                0.4294
## Genetic_DiversityLow                                                   0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.9722
##                                                                         
## (Intercept)                                                          *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.53 0.133 Inf      4.21      4.85
##  Medium              4.30 0.151 Inf      3.94      4.66
##  Low                 4.01 0.137 Inf      3.68      4.33
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.230 0.153 Inf 1.507   0.2874 
##  High - Low       0.523 0.140 Inf 3.746   0.0005 
##  Medium - Low     0.293 0.157 Inf 1.861   0.1503 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#Signifcant interaction 

3.4 Growth rate G6

head(data_phenotyping)
##   Population   Week  Block ID_Rep Divsersity Population.1 Box Initial.N   Start
## 1      High1 5-week Block3     52       High            1   1        30  7/8/21
## 2      High1 5-week Block3     53       High            1   2        30  7/8/21
## 3     High13 5-week Block4    129       High           13   1        30 7/15/21
## 4     High13 5-week Block4    130       High           13   2        30 7/15/21
## 5     High13 5-week Block4    131       High           13   3        30 7/15/21
## 6     High14 5-week Block4    132       High           14   1        30 7/15/21
##       End Larvae Pupae Adults    Lambda Genetic_Diversity Nb_ttx  p_adults
## 1 8/12/21      3     9     97 3.2333333              High    109 0.8899083
## 2 8/12/21      1     1      8 0.2666667              High     10 0.8000000
## 3 8/19/21      2     0     84 2.8000000              High     86 0.9767442
## 4 8/19/21      1     2    104 3.4666667              High    107 0.9719626
## 5 8/19/21      2     1     84 2.8000000              High     87 0.9655172
## 6 8/19/21      1     0     83 2.7666667              High     84 0.9880952
##      p_pupae    p_larvae He_remain  He_start    He_exp    He_end
## 1 0.08256881 0.027522936 0.9986008 0.9986008 0.9679591 0.9666047
## 2 0.10000000 0.100000000 0.9986008 0.9986008 0.9679591 0.9666047
## 3 0.00000000 0.023255814 0.9986008 0.9986008 0.9786327 0.9772634
## 4 0.01869159 0.009345794 0.9986008 0.9986008 0.9786327 0.9772634
## 5 0.01149425 0.022988506 0.9986008 0.9986008 0.9786327 0.9772634
## 6 0.00000000 0.011904762 0.9986008 0.9986008 0.9766045 0.9752381
#Without considering extinct populations 

#tapply(data_phenotyping$Lambda,data_phenotyping$Population,length)




#############################################################
################## Determine distribution ###################
#############################################################

## Poisson lognormal model
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
#summary(mlog)

# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
#summary(msqrt)

# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)


## Compare AIC
AIC(mlog, msqrt,mnorm)
##       df      AIC
## mlog   7 120.8153
## msqrt  7 156.0901
## mnorm  7 528.8983
#Square root a better fits to the data 
#But we compare different values: Growth rate and Growth rate+1

# ## Simulate data 
x.teo.log <- unlist(simulate(mlog))
x.teo.sqrt <- unlist(simulate(msqrt))
x.teo.norm <- unlist(simulate(mnorm))

# ## QQplot to compared Negative binomial and Poisson log normal distributions
qqplot(data_phenotyping$Lambda, x.teo.log,
       main="QQ-plot", xlab="Observed data", ylab="Simulated data",
       las=1, xlim = c(0, 10), ylim = c(0, 10), col='blue') ## QQ-plot
points(sort(data_phenotyping$Lambda), sort(x.teo.sqrt), pch=1, col='red')
points(sort(data_phenotyping$Lambda), sort(x.teo.norm), pch=1, col='green')
abline(0,1)
# ## Add legend
legend(0, 7, legend=c("Log", "Sqrt", "Normal"),
       col=c("blue", "red", "green"), 
       pch=1, bty="n")

# ## Normal distribution provides a good fit to the data


# QQ plot
# qqnorm(resid(mlog))
# qqline(resid(mlog))
# qqnorm(resid(msqrt))
# qqline(resid(msqrt))
# qqnorm(resid(mnorm))
# qqline(resid(mnorm)) #best fit

#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)

hist(residuals(msqrt),col="yellow",freq=F)  #Seems the better fit

hist(residuals(mnorm),col="yellow",freq=F) #Seems the better fit

#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
##     1-mle-norm 
## "not rejected"
## Compare AIC
AIC(mlog, msqrt, mnorm)
##       df      AIC
## mlog   7 120.8153
## msqrt  7 156.0901
## mnorm  7 528.8983
#SLog a better fits to the data 





#############################################################
##################        Analysis        ###################
#############################################################

mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
summary(mod0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##    Data: data_phenotyping
## 
## REML criterion at convergence: 514.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3747 -0.6242 -0.0778  0.5409  3.4451 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.2321   0.4817  
##  Residual               0.7480   0.8649  
## Number of obs: 186, groups:  Population, 57
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              2.569709   0.193717  13.265
## Genetic_DiversityMedium -0.873366   0.238587  -3.661
## Genetic_DiversityLow    -0.700898   0.222333  -3.152
## BlockBlock4              0.220484   0.214487   1.028
## BlockBlock5             -0.003695   0.249586  -0.015
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -0.460                     
## Gntc_DvrstL -0.587  0.363              
## BlockBlock4 -0.636  0.079  0.175       
## BlockBlock5 -0.592  0.089  0.232  0.451
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    5 532.53 548.66 -261.26   522.53                         
## mod0    7 520.90 543.48 -253.45   506.90 15.635  2  0.0004027 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE   df lower.CL upper.CL
##  High                2.64 0.135 42.2     2.31     2.98
##  Medium              1.77 0.198 51.5     1.28     2.26
##  Low                 1.94 0.177 58.2     1.51     2.38
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE   df t.ratio p.value
##  High - Medium    0.873 0.239 48.0  3.653  0.0018 
##  High - Low       0.701 0.223 52.5  3.143  0.0076 
##  Medium - Low    -0.172 0.261 53.0 -0.660  0.7873 
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
############ 
###### Predict
############
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_phenotyping$Genetic_Diversity),
                                   Block = "Block4")

data_predict_extinct$predict <- predict(mod0, 
                                        type = "response",
                                        re.form = NA,
                                        newdata = data_predict_extinct)

3.5 Remaining He difference

Rmisc::summarySE(data[data$Generation_Eggs==2,],
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd           se           ci
## 1              High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2            Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3               Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320
mfull <- lm(He_remain ~ Genetic_Diversity, data=data[data$Generation_Eggs==2,])
mless <- lm(He_remain ~ 1, data=data[data$Generation_Eggs==2,])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
##   Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
## 1     81 0.00631                                 
## 2     83 3.14572 -2   -3.1394 20162 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE df lower.CL upper.CL
##  High              0.9918 0.001698 81   0.9877   0.9960
##  Medium            0.6743 0.001840 81   0.6698   0.6788
##  Low               0.5404 0.001513 81   0.5367   0.5441
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.318 0.00250 81 126.829 <.0001 
##  High - Low       0.451 0.00227 81 198.470 <.0001 
##  Medium - Low     0.134 0.00238 81  56.200 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
Rmisc::summarySE(data[data$Generation_Eggs==2,],
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd           se           ci
## 1              High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2            Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3               Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320
#Calcul for end:
Rmisc::summarySE(data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",],
          measurevar = c("He_end"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N    He_end         sd          se          ci
## 1              High 27 0.9663114 0.01931641 0.003717444 0.007641317
## 2            Medium 18 0.6378048 0.03493729 0.008234799 0.017373906
## 3               Low 26 0.5090656 0.03443416 0.006753094 0.013908257
mfull <- lm(He_end ~ Genetic_Diversity, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
mless <- lm(He_end ~ 1, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_end ~ Genetic_Diversity
## Model 2: He_end ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     68 0.06009                                  
## 2     70 2.97522 -2   -2.9151 1649.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE df lower.CL upper.CL
##  High               0.966 0.00572 68    0.952    0.980
##  Medium             0.638 0.00701 68    0.621    0.655
##  Low                0.509 0.00583 68    0.495    0.523
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.329 0.00905 68 36.316  <.0001 
##  High - Low       0.457 0.00817 68 55.978  <.0001 
##  Medium - Low     0.129 0.00912 68 14.124  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.6 Growth rate and He

#Best model with genetic diversity
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)

mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod1 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod2 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod5 <- lme4::lmer(Lambda ~ exp(He_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
MuMIn::model.sel(mod0, mod1,mod2,mod3,mod4,mod5)
## Model selection table 
##       (Int) Blc Gnt_Dvr He_end log(He_end) exp(He_end)             family df
## mod1 0.8154   +          1.766                         gaussian(identity)  6
## mod5 0.3445   +                                 0.8311 gaussian(identity)  6
## mod3 2.5470   +                      1.259             gaussian(identity)  6
## mod0 2.5700   +       +                                gaussian(identity)  7
## mod2 2.0770   +                                        gaussian(identity)  5
## mod4 2.0770   +                                        gaussian(identity)  5
##        logLik  AICc delta weight
## mod1 -257.506 527.5  0.00  0.408
## mod5 -258.097 528.7  1.18  0.226
## mod3 -258.156 528.8  1.30  0.213
## mod0 -257.449 529.5  2.05  0.147
## mod2 -263.551 537.4  9.95  0.003
## mod4 -263.551 537.4  9.95  0.003
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Population'
#Makes sense: best model log He
x <- seq(1:100)
y <- seq(1001:1100)
plot(log(x),y)

#############################################################
##################        Analysis        ###################
#############################################################

mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_end) + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod4    5 532.53 548.66 -261.26   522.53                         
## mod3    6 522.65 542.01 -255.33   510.65 11.878  1   0.000568 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ log(He_end) + Block + (1 | Population)
##    Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## 
## REML criterion at convergence: 516.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3851 -0.6642 -0.0835  0.5336  3.4375 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.2567   0.5067  
##  Residual               0.7473   0.8644  
## Number of obs: 186, groups:  Population, 57
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 2.547160   0.201530  12.639
## log(He_end) 1.259386   0.351126   3.587
## BlockBlock4 0.229875   0.219317   1.048
## BlockBlock5 0.004275   0.254183   0.017
## 
## Correlation of Fixed Effects:
##             (Intr) lg(H_) BlckB4
## log(He_end)  0.655              
## BlockBlock4 -0.625 -0.154       
## BlockBlock5 -0.581 -0.197  0.446
#############################################################
##################        Predict         ###################
#############################################################

data_predict_lambda <- data.frame(He_end =
  seq(min(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end),
      max(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end), 
      0.01),
                                   Block = "Block4")

data_predict_lambda$lambda_predict <- predict(mod3, 
                                        type = "response",
                                        re.form = NA,
                                        newdata = data_predict_lambda)
 
  PLOT_He_expect <- PLOT_He_Final +   geom_line(data = data_predict_lambda, 
                                aes(x = He_end, y = lambda_predict),
                                linetype = "longdash", col = "grey40", size = 1.25) 


# # 
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"),   PLOT_He_expect,
#           base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)

3.7 Adults G6

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   1818.1   1834.3   -904.1   1808.1      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23572 -0.14532  0.02672  0.14825  0.41394 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.1739   0.4170  
##  Population        (Intercept) 0.1085   0.3294  
## Number of obs: 186, groups:  ID_Rep:Population, 186; Population, 57
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2943     0.0812  52.884  < 2e-16 ***
## Genetic_DiversityMedium  -0.5132     0.1438  -3.569 0.000358 ***
## Genetic_DiversityLow     -0.3047     0.1295  -2.352 0.018660 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.563       
## Gntc_DvrstL -0.625  0.356
# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population/ID_Rep),
#                 family = "poisson", data = data_5week)


#dispersion_lme4::glmer(m0) #dispersion ok 
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1    3 1826.3 1836.0 -910.15   1820.3                        
## m0    5 1818.1 1834.3 -904.07   1808.1 12.164  2   0.002284 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   1818.1   1834.3   -904.1   1808.1      181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23572 -0.14532  0.02672  0.14825  0.41394 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.1739   0.4170  
##  Population        (Intercept) 0.1085   0.3294  
## Number of obs: 186, groups:  ID_Rep:Population, 186; Population, 57
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2943     0.0812  52.884  < 2e-16 ***
## Genetic_DiversityMedium  -0.5132     0.1438  -3.569 0.000358 ***
## Genetic_DiversityLow     -0.3047     0.1295  -2.352 0.018660 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.563       
## Gntc_DvrstL -0.625  0.356
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean     SE  df asymp.LCL asymp.UCL
##  High                4.29 0.0812 Inf      4.10      4.49
##  Medium              3.78 0.1188 Inf      3.50      4.06
##  Low                 3.99 0.1011 Inf      3.75      4.23
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.513 0.144 Inf  3.569  0.0010 
##  High - Low       0.305 0.130 Inf  2.352  0.0489 
##  Medium - Low    -0.208 0.156 Inf -1.340  0.3728 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m1    4 1823.6 1836.5 -907.81   1815.6                       
## m0    6 1820.1 1839.5 -904.07   1808.1 7.4783  2    0.02377 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1    3 1844.6 1854.3 -919.30   1838.6                        
## m0    5 1838.5 1854.6 -914.26   1828.5 10.086  2   0.006455 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.8 Adults G2

m0 <- glm(Nb_adults ~ Genetic_Diversity + Block,  family = "poisson", data = data[data$Generation_Eggs==2,])
summary(m0)
## 
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson", 
##     data = data[data$Generation_Eggs == 2, ])
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.7057   -2.7621    0.2126    3.3983   11.4454  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              5.73890    0.01497 383.312  < 2e-16 ***
## Genetic_DiversityMedium -0.19408    0.01628 -11.920  < 2e-16 ***
## Genetic_DiversityLow    -0.19278    0.01460 -13.205  < 2e-16 ***
## BlockBlock4              0.10767    0.01579   6.817  9.3e-12 ***
## BlockBlock5              0.17743    0.01600  11.088  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2379.6  on 83  degrees of freedom
## Residual deviance: 2027.9  on 79  degrees of freedom
## AIC: 2667.2
## 
## Number of Fisher Scoring iterations: 4
m1<- glm(Nb_adults ~ Block,  family = "poisson", data = data[data$Generation_Eggs==2,])
anova(m0,m1,test="Chisq") # difference
## Analysis of Deviance Table
## 
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        79     2027.9                          
## 2        81     2241.4 -2  -213.52 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE  df asymp.LCL asymp.UCL
##  High               5.834 0.01051 Inf     5.809     5.859
##  Medium             5.640 0.01243 Inf     5.610     5.670
##  Low                5.641 0.01022 Inf     5.617     5.666
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate     SE  df z.ratio p.value
##  High - Medium   0.1941 0.0163 Inf 11.920  <.0001 
##  High - Low      0.1928 0.0146 Inf 13.205  <.0001 
##  Medium - Low   -0.0013 0.0161 Inf -0.081  0.9964 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates